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Data Assimilation for Convective Cells Tracking on 


Meteorological Image Sequences 

Claire Thomas* Thomas Corpctti 'and Etienne McmiiV 

Abstract 

This paper focuses on the tracking and analysis of convective clouds systems from Meteosat Second 
Generation images. The highly deformable nature of convective clouds, the complexity of the physic pro¬ 
cesses involved, but also the partially hidden measurements available from image data make difficult a direct 
use of conventional image analysis techniques for tasks of detection, tracking and characterization. In this 
paper, we face these issues using variational data assimilation tools. Such techniques enable to perform the 
estimation of an unknown state function according to a given dynamical model and to noisy and incomplete 
measurements. The system state we are setting in this study for the clouds representation is composed of 
two nested curves corresponding to the exterior frontiers of the clouds and to the interior coldest parts (core) 
of the convective clouds. Since no reliable simple dynamical model exists for such phenomena at the image 
grid scale, the dynamics on which we are relying has been directly defined from image based motion mea¬ 
surements and takes into account an uncertainty modeling of the curves dynamics along time. In addition 
to this assimilation technique, we show in appendix how each cell of the recovered clouds system can be 
labeled and associated to characteristic parameters (birth or death time, mean temperature, velocity, growth, 
etc.) of great interest for meteorologists. 
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1 Introduction 

Convective cells are atmospheric events that are known to be associated with hazardous consequences, such 
as strong wind drafts, lightnings, heavy rainfalls, hails or even tornadoes. Over tropical areas such as central 
Africa, convective cells produce most of the rain during the monsoon period. They are also indirectly linked to 
droughts and floods, which might afflict this area. Their analysis and forecasting are thus of the utmost interest 
for meteorologists and forecasters. 

The measurement of the convective cells parameters can be done either by the use of conventional probes 
or through satellite data. However, over remote areas (like the Amazonian forest for instance), conventional 
sensors such as radiosonde, rain gauges, radars or lightning detectors have a limited coverage. Satellite in¬ 
formation constitutes therefore an appealing alternative for the study of convective atmospheric activity over 
these particular regions. 

Except for large convective systems, the study of the convective cells development is particularly intricate 
since they are generally very sudden, short-lived and constitute highly deformable objects. As a consequence, 
the study of the convective activity requires high spatial and temporal resolution for precise diagnosis and fore¬ 
casts. The new generation of geostationary space-borne sensors provides a wider variety of image channels 
with narrower bandwidths and an increased spatial and temporal resolution suitable for the analysis of synop¬ 
tic to mesoscale phenomena. In this study we will rely on data of the Infrared (IR) and Water Vapor (WV) 
channels of the Meteosat Second Generation satellite (MSG-2). Compared to the first version, this sensor has 
been improved both on spatial resolution (1 pixel corresponding to 3km instead of 5km in the first generation) 
and temporal sampling rate (15min between two images instead of 30min in the previous system), in response 
to users requirements. 
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As introduced previously, this paper aims at studying methods for the reconstruction of a convective clouds 
system trajectory. Despite of the spatial and temporal resolution of satellite images, this task remains challeng¬ 
ing, as images constitute only a 2D discrete sequence along time of a 3D phenomenon. Besides, as convective 
cloud systems exhibit a complex dynamics in which cells may split or merge, the 2D shape of such clouds 
system may change drastically between two successive frame instants. The reconstruction of a continuous 
trajectory to infer the whole history of a convective system in term of cells birth/death time, or of their activity 
phases constitutes hence an intricate issue that can hardly be handled by running “decoupled” segmentation 
processes on successive images. 

In the past, several attempts have been proposed in the meteorological community for an automatic detec¬ 
tion and tracking of these phenomena from satellite image data [1,2,17,29-31,33,34,39]. The first detection 
approaches were based on temperature threshold of IR images. These simple techniques permit to detect the 
cold core of convective cells. They, however, do not enable a reliable extraction of the cells contours, which 
are very sensitive to the threshold value. A prediction of the cells location roughly identified to the cold blob 
gravity centers has been then carried out through correlation techniques. Assuming an effective overlapping 
between the predicted cells and the detected ones at the currant instant, such techniques provide successive es¬ 
timations of the cells core average location. Velocity and sometimes others parameters (such as the ellipticity 
factor or the distribution of the temperature gradient within the cells cores) are also provided. 

In computer vision, a primal strategy to delineate the contours of complex objects relies on the implemen¬ 
tation of partial differential equations encoding the evolution of parametric or non parametric curves toward 
the minimizers of an energy functional [8, 9, 24], The corresponding functionals include generally a data 
term representative of a photometric distribution characterizing the object of interest and some regularity con¬ 
straints on the curves. These non linear deterministic minimization strategies are in practice very sensitive 
to the initial curve location and topology. Among the different extensions of those methods, the Eulerian 
“level-set” approach [35,43] has been specifically proposed to alleviate such shortcomings. In this framework, 
the contour’s shape is represented as the zero level-set of an implicit higher-dimensional scalar function. The 
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evolution of this implicit surface describes the contour’s evolution and naturally enables the handling of topol¬ 
ogy changes. Such techniques have been already proposed in a meteorological context [10,39,51]. However 
they generally remain static detection/segmentation approaches applied (quasi-)independently on the images 
of the sequence with no temporal consistency guaranty on the recovered shape sequences. Time consistency 
is brought considering the estimation of the whole curves trajectory on the basis of a system of equations 
that couples an evolution model representative of the target dynamics and the complete image sequence. This 
constitutes a data assimilation problem that can be handled either relying on stochastic filtering or variational 
assimilation techniques. Stochastic filtering techniques implemented through particle filters have been pro¬ 
posed for the tracking of parametric curves [7] or level set [40]. These techniques, which have the advantage 
of being recursive are nevertheless limited by difficult sampling issues in high dimensional state spaces. This 
limitation refered in the litterature as the curse of dimensionality, requires in practice the use of simple para¬ 
metric descriptors for the curve (which prevents any topological changes) or the definition of simplified linear 
motion models. This latter simplification makes impossible the reconstitution of complex deformations and 
can hardly describe accurately the dynamics of convective cells. Compared to stochastic filtering techniques, 
the variational data assimilation framework [27], provides a convenient procedure to handle high dimensional 
data assimilation problems. This technique belongs to the class of smoothing approaches, that operates on the 
whole sequence through forward and backward integrations of two coupled systems of pde’s (the dynamical 
model and its adjoint). This process supplies a single trajectory, solution of the specified dynamics that fits at 
best a sequence of observations related to the state variable of interest. 

In this work, we aim at specifying an image based variational data assimilation procedure [11,36-38] for 
the tracking of a convective clouds system. The definition of such a scheme requires to establish an observation 
operator linking the unknown curve with the available data and to define a dynamics that i) represents, as 
accurately as possible, the phenomenon of interest and ii) allows an efficient estimation of the unknown state 
variables. In our case this specification should describe the evolution along time of the 2D projection of the 
external frontiers of a convective cell and of its core. An accurate dynamics defined at the image grid scale that 
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accounts for the different physical processes involved in those meteorological phenomenon is far to complex 
to be included in a practical image based assimilation process. A simplified “imperfect” model has then to be 
considered. Such assimilation scheme referred as weak constrained variational assimilation relies generally 
on a physical model simplification defined up to an uncertainty variable. In this work, we will comply with 
the same principle. However instead of sticking to a simplification of the cells dynamics, which would be in 
its own a difficult issue (as to our knowledge no satisfying evolution models are for the moment available), 
we will rely on a generic dynamics driven from image based motion measurements. This generic model will 
allow us to describe the evolution of two nested curves representing the cell cores and their external contours. 
The simplicity of this model will be paid in return by the fact that as we do not have any dynamical models for 
the flow within or in vicinity of the cells, it will be impossible to set up a denoising data assimilation process 
of the velocity measurement as proposed in [11,36]. Nevertheless, since the model considered incorporates an 
uncertainty related to the curve dynamics, we will show that the proposed model is robust to inaccurate motion 
measurements and does not require any specific assimilation of the flow velocity fields. 

The structure of the paper is organized as follows: section 2 briefly describes the data used and justifies the 
choice of specific MSG channels; section 3 outlines the main principles of variational assimilation. In section 
4, we describe the proposed data assimilation procedure for the tracking of two nested level sets delineating 
respectively the core and the boundary of the convective cell. The section 5 presents some experimental results 
and discusses the tracking robustness. Section 6 draws some conclusions on this work. Finally, the appendix 
opens new perspectives on the use of the estimated level set for the temporal analysis of convective cells. 

2 Dataset 

2.1 Channel for the convective cell detection 

The first issue consists in describing properly the clouds designated as “convective cell” in order to define an 
adequate geometrical description together with photometric signatures that characterize them. In the literature 
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the term “convective clouds” carries some ambiguities as it is employed to designate a multitude of very dif¬ 
ferent meteorological phenomena such as: derecho, thunderstorm, multicell complex, supercell, squall line or 
convective cell. As explained by Morel and Senezi in [33], convective cells, or Mesoscale Convective Systems 
(MCS), are convective phenomena that take place at meso-o scales (250-2500 km) or upper meso-/f scales 
(25-250 km). They have been studied for decades. In our study, we focus on the meso-/3 cells, and rely on 
the spatial description proposed for those cells by Houze in [21]. A convective cell is a cloud producing low 
temperatures, and composed of two distinct regions: the actual convective part which consists of intense, cold¬ 
est, vertically extended cores, and the stratiform region characterized by a more uniform texture and lighter 
precipitation. The stratiform area is partly produced by the dissipation of older convective cells and partly 
produced by broader sloping mesoscale layer ascent. The horizontal pattern formed by the convective core 
and the stratiform cloud part exhibit large variations along time and constitute both important indicators of the 
internal cloud dynamics. In this work, we will therefore describe a convective cell as two nested curves in 
order to distinguish the stratiform area from the actual convective part of the cloud. 

The discrimination of convective cells is usually performed using a temperature threshold on IR imagery, 
since those events constitute the coldest events observed by satellite. In order to isolate these cloud systems, 
the image luminance needs to be transformed into temperatures. In general, physical temperatures of the sky 
can not be directly inferred from satellite imagery [45]. However relying on a black body assumption an 
equivalent blackbody brightness temperature can be extracted for top of clouds observed through IR channels. 
This assumption is valid only for thick clouds like the convective cells. However, other elements such as gases 
and aerosols located between the cloud of interest and the satellite sensor may absorb or diffuse a portion of 
the radiation emitted by the cloud. As a consequence, due to strong water vapor absorption in the WV 6.2 
and WV 7.3 channels, these channels cannot be used for temperature determination, as opposed to IR bands, 
which are thus commonly adopted for the detection of convective system. 

Relying on previous studies [2,21,31,33,45], we have selected the IR 10.8 channel for brightness temper- 
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ature estimation. Two thresholds (Ti = —20 °C = 253A' and T 2 = — 60"C = 2 13 A) are used to separate 
respectively the convective cells and the core from the other clouds. 

Let us recall that, following the formation process image, the brightness temperatures are computed from the 
IR luminance values using: 


T b = 


C 2 v c 


(2# + l) 


— B 


Ain 


( 1 ) 


where C\ and C 2 are constants respectively equal to 1.1910410 -5 mfLrn _2 sr _1 cm -4 and 1.43877A” cm. 
Parameters A and B differ for each satellite/channel and their values are available on the EUMETSAT web¬ 
site http : // www.eumetsat.int for infrared and water vapor channels (IR and WV). The radiance R (in 
mW m~ 2 sr -1 cm) is received by the detector and is linearly related to the luminance value I by: 


R — Rq T n /. 


( 2 ) 


where a and Rq are constants that might slightly change from day to day and whose values are provided by 
EUMETSAT. Let us note however that some very localized clouds associated to shorter lifetime (between 30 
minutes and one hour) and a weaker convection activity exhibit brightness temperature below the threshold 
of —20°C. Even if these clouds reach also the tropopause, they do not belong to the class of Mesoscale 
Convective Systems. To discard those isolated cloud storms, it is possible to discriminate them based on 
the lightnings occurences [33]. However, this requires additional data. In our study, this discremination will 
be performed naturally through the proposed assimilation process. Let us now present the data used for the 
motion estimation technique. 


2.2 Channel for the optic flow estimation 

As further presented in section 4, the assimilation process we built requires the estimation of velocity fields 
describing the top of the convective motion. On the basis of comparisons studies leaded on atmospherical wind 
field estimation from meteorological satellite images [16,26,45], we selected the WV channels to proceed the 
wind estimation. Let us point out that an important characteristic of pure water vapor structures is that they 
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represent a mixture of radiation originating from different altitudes, located between pressure levels 700 and 
200 hpa, (i.e. altitudes between 3 and 12 km). Consequently, WV channel provides on pure water vapor 
area an integration of structures at different heights. In the case of the convective cells, thunderstorms have a 
thickness of several kilometers and their wind motion information might strongly differ at the cloud’s top and 
at its base. The use of WV channels for the determination of an average wind field related to the convective 
cells is hence pertinent. The wind estimation may however not be relevant in areas without clouds and with a 
poor water vapor content [26]. This has to be taken into account in the assimilation process. 

3 Data Assimilation 

In this section we will describe the general principles of the assimilation scheme we devised in this study. 

3.1 Definition 

Variational data assimilation is a technique derived from optimal control theory [28] to recover a state vari¬ 
able’s trajectory from a sequence of measurements. Opposite to sequential Bayesian filters, which share the 
same aim, this framework allows to handle high dimensional systems and is thus intensively used in environ¬ 
mental sciences [6,27,47]. In this study, we rely on this framework to estimate the complete trajectories along 
an image sequence of closed curves representing convective clouds systems. 

The data assimilation setup is basically composed of an ideal dynamics of the system variables, an initial¬ 
ization of the system’s state and an observation equation that relates the system variables to some measure¬ 
ments: 

d t X + M(X) = u(x, t) 

' X(x,t a ) = X 0 (x) + r/(x) (3) 

Y(x,t) = H(X (x, t)) +e(x,t). 

The right hand side of the first equation describes, through a differential operator M, the evolution of the state 
function X(x,t) defined over the image plane 12 and the whole image sequence time range [t, 0 ;tf\. In our 



case, the components of this function will be composed of two implicit scalar surface functions representing 
two nested closed curves. This will be further detailed in the next section, but for sake of generality, let us first 
consider a general state function. The second equation sets up an initial condition for the state function through 
a given initialization X a (x). The last equation links an observation function Y(x, t) to the state function. This 
relation is formalized through a differential operator EL The observation function is usually composed of noisy 
and sparse measurements. In these three equations v , // and e are time varying Gaussian noise functions defined 
on the whole image plane. Since in our case the observation depends directly on the luminance function, 

I(x,t), whereas the dynamics depends on motion measurements estimated from the luminance spatial and 
temporal gradients, (V/, dt.I), we assume that v and e are conditionally independent with respect to the image 
sequence. Besides, in order to simplify the estimation scheme and in order to fit to simplified variational 
assimilation schemes we will assume these noises are correlated in space but uncorrelated in time. This 
restriction could easily be dropped in the following but at the price of a much higher computational cost. Noises 
v, tj and e are thus respectively associated to definite positive endomorphisms Q(x, x')6(t — t 1 ), B(x, x') and 
R(x,x')5(t — t') refered as covariance tensors. These tensors represent the errors involved in the different 
components of the system (i.e dynamical model errors, initialization errors and measurement errors). 

3.2 Penalty function 

The system of equations (3) could as well be specified through three conditional probability densities p(X(., t)|X(., t 0 ), /), 
p(X(.,t 0 )\X 0 (.)) and p(Y (., t)|X(., t), I) related respectively to 

i) the distribution of the state trajectory X(., t) given the initial state value (., t 0 ) and the whole image 
sequence/ = {I(.,t),t G 

ii) the distribution of initial state value given the initial condition X a (.) and the image sequence, and 

Hi) the distribution of Y(., t), the complete set of measurements, given the state X and the image sequence. 
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For the system we are considering the second distribution is Gaussian whereas the first and the third ones 
are Gaussians only if their associated operator (observation and dynamic) are linear. In practice the dynam¬ 
ical operator is highly non linear and due to for instance 2D-3D projections or indirect measurement of the 
quantity of interest, the associated observation operators we have to cope with in image are often non-linear. 
Generally, it is considered that the less linear these operator, the more their distributions depart from Gaussian 
distributions. 

As in any stochastic filtering problem, we aim here at estimating the conditional expectation of the state 
trajectory given the whole set of available measurements and the complete image sequence. A maximum a 
posteriori estimation consists in estimating the mode of the distribution p(X\Y, X 0 .1). Assuming Gaussian 
approximations of the pdf involved and uncorrelated errors in time, this leads to the minimization of a quadratic 
functional defined as: 

i rtf t i As 

J(X,u, V ) = -j^ \\Y-W{X)\\%- 1 dt+-\\X{t 0 )-X 0 \\%- t +-J t ||—+M(X)|| l^dt, (4) 

where we have introduced the norm ||JA||associated to the scalar product in L 2 (fl): 

(X, A _1 Y) — f j X(x)A(x,x')Y(x')dxdx'. 

Jq Jq 

When the assimilation relies on an exact version of the dynamics, the last term of this functional (the 
model part) disappears. In that case, the functional then depends only on the initial condition, and comes to 
an initial value control problem [27,47]. This scheme is commonly referred in the literature as variational 
data assimilation and abbreviated as “4Dvar” whereas associated to dynamics with uncertainty variables, it is 
refered as “weak-constraint” variational data assimilation methodology [4,5,15,48,49], Both schemes rely 
on the same adjoint optimization principles [27,47] to minimize the functional. In the following we give a 
presentation of this principle. 

A minimizer X of functional J is also a minimizer of a cost function J(X + f30(x,t)), where 9(x,t )) 
belongs to a space of admissible functions and [3 is a positive parameter. In other words, X must cancel out 
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the directional derivative: 


d x J{6) = lim 

/ 3—0 


dJ{X + f39(x, t)) 
d/3 


= 0 . 


The computation of this gradient can be realized using an “adjoint formulation”, as explained in the next 
section. 


3.3 Functional gradient 


Introducing an “adjoint variable” A 


X(x,t)= [ Q- 1 (d t X + M{X))dx', 

Jn 


(5) 


as well as the Gateaux derivatives at point X of the operators M and H (called the linear tangent operators ), 
denoted and d x H respectively and defined by: 


.. dM(X + 00) _ 

lim- = d x M6, 


/3—>o d/3 

the directional derivative of functional J ( X , v, rf) reads: 


(6) 


Km ^ = (0( to ), B~ 1 (X(x', t 0 ) - X o) ) - J ( d x m , R-\Y - H(X))) dt + J (d t 0 + d x M6 , A) dt (7) 

Considering the three following integrations by parts, one can get rid of the partial derivatives of the admissible 
function 9 in expression (7), i.e. 


I(d t 9, A) dt = (0(tf),\(tf))-(O(t o ),\(x,t o ))- / ( 9,d t X)dt , 


{d x MS, A) dt = / (9, d x M*X) dt, 


> x m,R- 1 (y-m{x)))dt= / (9,d x m*R- 1 (y -u{x)))dt. 


( 8 ) 


(9) 


( 10 ) 
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In the two last equations, we have introduced adjoint operators cf^M* and dxB* as compact notations of 


the integration by parts of the associated linear tangent operators. Gathering all the elements we have so far, 
equation (7) can be rewritten as: 

dxJ{0) = + {6(t 0 ),B~ l {X{t 0 ) - X 0 ) - A (t 0 ))dt+ 

,(-d t \ + dxM*\)-d x B.*R- 1 (Y-W(X)))dt. (11) 

This functional gradient must cancel out to satisfy to optimality conditions. This yields a “forward/backward" 
integration of the linear tangent dynamical system and its adjoint. 

3.4 Forward/backward equations 

Since the functional derivative must be null for any arbitrary independent admissible functions, the right mem¬ 
bers of the three scalar products of expression (11) must be identically null. It follows a coupled system of 
forward and backward PDE’s with initial and final conditions: 


A(*,t / ) = 0, (12) 

X(x,t 0 ) = [ B~ 1 (x,x')(X(x',t 0 ) - X 0 {x'))dx', (13) 

Jn 

dtX(x,t)+M.(X(x,t))= f Q{x,x')X(x',t)dx\ (14) 

J n 

d t X-d x M*X=- [ d x m*R- 1 {x,x')(Y -M{X))dx'. (15) 

Jn 


The forward equation (14) corresponds to the definition of the adjoint variable (5) and has been obtained 
introducing Q, the pseudo-inverse of Q ~ 1 , defined as [6]: 

f Q(x,x')Q~ 1 (x', x”)dx' = 8(x — x")8(t — t"). 

Jn 

Let us remark that for Q = 0 we retrieve the case of a perfect dynamics. Otherwise, equation (12) constitutes 
an explicit final condition for the adjoint evolution model equation (15). This adjoint evolution model can be 
integrated backward from this final condition assuming the knowledge of an initial guess for X to compute 
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the discrepancy Y — H(X). To perform this integration, an expression of the adjoint evolution operator is 
required. Let us recall that this operator is defined from an integration by part of the linear tangent operator 
associated to the evolution law operator. The expression of such an operator can be obtained analytically or 
from the discrete scheme associated to the linear tangent operator. As a matter of fact, the adjoint of the linear 
tangent operator discretized as a matrix simply consists in the transpose of that matrix. 

Knowing a first solution for the adjoint variable, an initial update condition for the state variable can be 
obtained from (13) and a pseudo inverse expression of the covariance tensor B. From this initial condition, 
equation (14) can be finally integrated forward. 

It should be noted that equation (14) provides at convergence an expression of the model error: 

v(x,t) = / Q{x, x')\{x', t)dx'. (16) 

Jn 

The knowledge of the adjoint variable enables thus to estimate the error associated to the state variable evo¬ 
lution law. This may be sometimes useful to validate or invalidate a tracking result, or at least to infer some 
confidence measures on the result obtained. 

This scheme somewhat differs from traditional weak constraint data assimilation. As a matter of fact, 
in geophysical applications, for computational reasons it is generally impossible to consider an uncertainty 
variable of the same dimension as the state variable. To control this variable a dimension reduction is generally 
applied or an explicit deterministic evolution model is imposed to the uncertainty model in order to come back 
to an initial value control problem [49], In our case, the state space is much smaller than in those applications 
and a full uncertainty model can be applied. 

3.5 Incremental state function 

The previous system relies on a Gaussian assumption ofp(X(.,t)\X(., t 0 )) and p((Y., t)\X(.,t)) (and implic¬ 
itly thus on linear assumption of the model and observation operators). Besides, as due to the non linearity of 
the operators involved in the minimization described previously, this procedure constitutes a direct non-linear 
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least squares estimation of the state variable. This process is likely to converge toward a bad local minimum 
if the initial guess is too far from the sought solution. The previous minimization can be slightly modified 
to follow a Gauss-Newton incremental methodology [14,48], which relies on successive linearization of the 
operators (and therefore better fits to the Gaussianity assumption). This scheme is built introducing a function 
of state increments linking the state variable and a background solution: SX = X — X Q . The linearization of 
the model operator M around the background X q: 

M(X) = M(X 0 ) + d Xo M(fiX), 

enables to split equation (14) into two pde’s with an explicit initial condition performing a model integration 
from a first guess solution X 0 (x)\ 


X(x,t a ) = X 0 (x), 

(17) 

d t Xo + M(A' 0 ) = 0, 

(18) 

d t SX + dx 0 M5X = f Q(x, x')\{x', t)dx'. 

J n 

(19) 


Combining equations (12-13) and (17-19) leads to the final assimilation algorithm. The method first consists 
of a forward integration of the background X a with a perfect system dynamics (18). The current solution 
is then corrected by performing a backward integration (12, 15) of the adjoint variable. The evolution of 
A is guided by a discrepancy measure between the observation and the estimate: Y — M(X ). The initial 
condition is then updated through equation (13) and a forward integration of the increment SX is realized 
through the equation (19). The overall process is iteratively repeated until convergence. A sketch of the whole 
process is summarized in figure (1). As a Gauss-Newton optimization process this non-linear minimization 
may fail to converge if the initial condition of the incremental function is too far from the sought solution. In 
our case, we will show experimentally that even for an incorrect initial trajectory the system still show good 
convergence property. This nice behavior is due to the observations we considered, and that enables to correct 
efficiently the initial trajectory. These observations are spatially highly resolved and allow getting noisy but 
good representations of the clouds system at discrete instants. The assimilation process manage to removes the 
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Let X Q (to) — X a . 

(i) From X 0 (t 0 ), compute X a (t), Vi G 
]t 0 ,tf[ with a forward integration of 
system (18). 

(ii) X Q (t) being given, realize a back¬ 
ward integration of the adjoint vari¬ 
able with the system (15). 

(iii) Compute the initial value of the incre¬ 
mental function (13). 

(iv) From SX(t 0 ), compute SX(t), Vi G 
]t 0 ,tf[ with a forward integration of 
system (19). 

(v) Update X = X + SX. 

(vi) Return to (ii) and repeat until conver¬ 
gence. 


Figure 1: Incremental variationnal data assimilation process 

noise or the false alarms associated to these observations and provides a continuous trajectory of the convective 
clouds system. 
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4 Application to convective cells 


This section presents the principles of the variational data assimilation scheme we propose for the tracking of 
a convective cells system represented by a set of nested curves. In is organized as follows: 

• the section 4.1 introduces the system state to assimilate; 

• the section 4.2 is devoted to the dynamical model; 

• the section 4.3 deals with the initial condition; 

• the section 4.4 presents the observations and 

• the section 4.5 is concerned with the definition of the error covariance matrices. 

We now turn to the definition of the system state. 

4.1 System state: convective cells representation 

Among the different possibilities available to represent closed curves, the level sets formalism [35,42] has 
proves its ability to manage in an elegant way the changes of topology (split, merge, ...) that are likely to 
appear. As the convective cells often undergo rapid changes in their topology and shape, such a formalism is 

well adapted to our problem. Within such framework, the curve at time t, Ct{p) = C(p, t) : [a, b] x RH-» R 2 , 

is implicitly described by the zero level set of a scalar function cf>(x, t) : O x RH -> R: 

C t {.) = C(.,t) = {s G fl | <p(x,t) = 0}, 

where stands for the image spatial domain. The surface is driven by the contours dynamics and is chosen so 
as to have for instance positive values inside the curve and negative values outside. A common choice for the 
implicit function is the signed distance function but any other surfaces whose level sets fits the curves of interest 
is possible. This representation has the great advantage to allow describing through a single implicit surface a 


16 



set of non intersecting closed curves. Besides, geometric properties of the curve can be directly obtained from 
the implicit surface. The inward unit normal and the mean curvature are for instance respectively given by: 


V0 

~W\ 


and 


k = div 



( 20 ) 


In order to encode the description of a convective cloud system encapsulating a set of convective cells, the 
system state X we used is then composed of two implicit surfaces: 


• The first one, (if, is devoted to the description of the external delimitations of the convective cells, C e (t). 
These curves will be associated to a first level of brightness temperature threshold T e . 


• The second one, (f> c represents the cold core region of each cells, C c (t). It is related to a colder threshold 
value T c of the brightness temperature. 


The first level set is defined on the whole image domain O: 


(jf{x, t) : fl x Rd-> R, 


whereas the second one is defined inside the domain encapsulted by the convective cells contours, X = 
{a; | <p e (x, t) > 0} (i.e. the domain delineated by the zero level set of <p e ): 

4> c (x , t) : X x RH-» R. 


The contours of the convective cores, represented by the zero level set of the second implicit surface, are 
thus constrained to lie inside the domain delineated by the external frontiers of the convective cells set. Both 
implicit surfaces constitute the components of the state variable of our data assimilation problem. Let us now 
define the dynamical model. 


4.2 Dynamical model 

In this section, we first introduce (in 4.2.1) a common way to represent the evolution of implicit surfaces, 
where these latter are assumed to be driven by a smooth external motion field without any uncertainties. As 
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will be shown, this assumption needs to be relaxed and we then propose a new and original dynamical model 
for the two nested curves where the motion uncertainties are explicitly taken into account (in 4.2.2). The 
resulting adjoint and linear tangent operators, some associated numerical aspects and the way to obtain the 
external motion fields are successively presented in sections 4.2.3, 4.2.4 and 4.2.5. 

4.2.1 Implicit surface evolution without motion uncertainties 

As commonly done when one deals with implicit representations, we assume that the curves of interest Q, i = 
{e, c} are driven by a smooth external motion field w(x. t ) along the normal n* and diffuse according to mean 
curvature flow [35,42]: 

nc 1 

-—t = (w l ■ n l + an l )n l , i = {e, c} (21) 

at 

where a and k respectively represent a diffusion coeficient and the curvature of the closed curve. As (f> 1 (Q ( p ), t ) 
0, the implicit surface dynamics reads: 

8C i 

d t <t> i +V<l> t -^= 0. (22) 

at 

Replacing in this equation the expression of the curve velocity (21) and of the normal (20), and extending this 
evolution law to the whole plane, we obtain finally the following transport equation: 

d t <l> i + V<t> i -w i + <TK\V<f> i \=0, (23) 

where we have considered a velocity field w l (x,t ) defined as an extension to the whole plane of the con¬ 
tour C‘(.,t) velocity. This representation provides an Eulerian representation of the evolving contours. As 
such, it allows to avoid the ad-hoc regridding processes of the different control points required in any explicit 
Lagrangian - spline based - curve description. 

However, one can remark that this dynamics relies on (i) an ideal noise free external motion field and 
requires (if) a smoothing parameter on the curves to guarantee a viscosity solution [42], This parameter is not 
easy to fix, as it is not related to any physical quantity of the object of interest. We will see in the following 
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section how this dynamics can be adapted to include some uncertainties on the motion field and how this 


provides us as a by product a more meaningful expression of the curve smoothing term. 


4.2.2 Implicit surface evolution with motion uncertainties 


In this section we introduce uncertainties on the motion field w of equation (21) that drives the two curves. To 
that end, we will rely on a stochastic formalism. For an introduction to stochastic calculus, we refer the reader 
to [23], 

We assume that the curves C\ are two stochastic processes driven by an Ito diffusion: 

dC\ — (w 1 ■ n l )n l dt + a n n l dB " + a T n l± dBl . (24) 


Here ri l dB n and n l± dB T are two independent Brownian motions directed along and perpendicularly to the 
curves normals and cr n , a T are the uncertainties respectively related to the normal and the tangent of the curve. 
In this model, the curve is transported by a deterministic drift associated to the velocity field w (first term 
of relation (24)) that is mitigated with isotropic Gaussian incertitudes -along the normal and tangent- whose 
covariances grow linearly in time (second and third terms of (24)). As previously, the curve location is defined 
from the evolution of an implicit surface cf> l (x,t). As this surface depends on the stochastic process C t , its 
differential must be calculated using stochastic calculus differentiation rules, the so called Ito formula [23], 
We get the following stochastic partial differential equation (in the following, for sake of clarity, we will drop 
the curve’s index i = {e, c} unless explicitly needed): 


( d t q ) + w • V0 + - 




E 

}={a 


,y} 


d 2 (j) 

dxidx. 


■d{C?\Ct j ))dt + a n \\7(j)\dB? = 0, 


(25) 


where Cf and C‘( correspond to the coordinates of the curve-points in the Euclidian plan at time t. As pre¬ 
viously, we have introduced velocity drifts w(x) and Brownian motions B™(x) which are extensions on the 
whole plane of the curves drift and noise. For the drifts, we considered a unique smooth velocity field defined 
on the whole plane for both curves. This velocity field will be further detailed subsequently in section 4.2.5. 
As for the noise, we considered a noise function defined from a set of independent Brownian noise variables 
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localized on the points of a fixed discrete grid and a set of gaussian smoothing functions, ip v , with variance v 
and centered on the same points: 

B t {x) =^2<p v (x- Xi)B?(xi). (26) 

i 

For a number of point tending to infinity, this function admits a limiting covariance tensor defined as: 

C(x,y,t,s) = ip 2 u{x-y)S(t-s). (27) 


The notation 


(X t ,X t ) = lim Ati ^((X ti+1 - X u ) T {X ti+1 - X u )) 


u<t 


involved in (eq. 25) denotes the quadratic variation of X t . Following stochastic calculus properties, this 
quantity is null for continuously differentiable functions and the following relation holds: 


(J <J t dB t , J a t dB t ) = J a 2 dt. 


The term d{X t ,X t ) = dXjdX t is computed according to the rule dtdB t = dtdt = 0 and dB t dB t = dt. In 
our case, we get: 


^ 

d{C!,C?) = 


d{c;.Cf) = 


|V0| 2 

On - 
|V0| 2 


dt. 


(28) 

(29) 

(30) 


Introducing the surface normal expression, the diffusion driving the implicit surface evolution reads finally: 


(d t (j) + w • V0 ■ 


2|V0| 


2 O’xzOl'/’z + ^2 4>y) + <t>yy{vi<t>v + <?2<l>x) + 2 On - a T )<t> x <l>y<f>xy))dt 


<r n \V<j)\dB = 0. (31) 


With the expression of the mean curvature: 


k = curv(^) 


b4 ( ^- v ^v^v«. 
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we get, 

(5*0 + M(0))df + cr„|V0|dB = 0, (32) 

with 

2 2 

M(0) = u> • V0 + ^/t| V0| + ^_V0 T V 2 0 V0, (33) 

where V 2 0 denotes the Hessian matrix. It can be observed from (33) that if both incertitudes have the same 

strength (i.e. a n = a T ) this model takes a particular simple form, which reads: 

2 

(dt4> + u> ■ V0 + -^-KA<f>)dt + <j n \V(f>\dB = 0, (34) 

The dynamical model (32-33) constitutes a general stochastic process allowing to guide a curve through 
an implicit surface. Nevertheless, as in our work, we wish to stay within a deterministic framework, we will 
assume a Gaussian approximation of the state variable probability distribution. As the expectation of the 
dynamical model is null we get the following approximation for the model pdf: 

p(0t|0o) otexp- ^ || d t (j> + M(0)||q_i (35) 

where from (27) the covariance tensor is 

Q(x,y,t,s) = al\V<j)(x,t)\ x |V0(y, s)\C(x, y)S(t - s). (36) 

This prior defines our complete evolution model. Let us note that the inverse of this covariance tensor is not 
required in the assimilation process of relation (19) and does not have to be computed. Compared to the model 
in (23) where no motion uncertainties were taken into account, the new formulation in (32-33) is more generic 
since it includes both uncertainties along the normal and the tangent of the curve. It is interesting to observe 
that on the deterministic part of the model (first term of relation (32)), the motion incertitude along the curves 
tangent corresponds to a smoothing of the curve and leads to the mean curvature flow introduced in (23). The 
uncertainty along the curves normal generates also a smoothing on the basis of the Hessian V 2 0 of the levelset. 

In practice, the constant a n and o> are fixed from an external observed displacements field between two 
consecutive images. More precisely, this field d(x,t) : I 2 x I -t R 2 , given in practice by a procedure 
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described in section 4.2.5 is assumed to be constituted by noisy measurements on the normal and the tangent 
of the photometric level-lines and of the form: 

d(x, k) = w(x, k) + a n (B(x, k + 1) - B(x, + aT ^^ x ’ k + l)~- B(x, k)) |yj_j| • (37) 

We are therefore making here the assumption that the observed velocity noises and the curve noises are 
collinear and have the same variances. Furthermore, we assume that the smooth drift velocity field, used 
to drive the curves, is such that w(x, t) = E(d(x, k))/dt and we approach this expectation through a local 
spatial mean over a small neighborhood V (x) of point x: 

w{x,k) = 1 ( d(x u k )). (38) 

1 ^ Xi&V(x) 

The variances of this observed displacement projected along the normalized curve normal and tangent provide 
an estimation of the noise variances <j„ and n't. Let us now turn to the constitution of the tangent linear 
operator and of its adjoint necessary to the assimilation process application. 


4.2.3 Adjoint and linear tangent operator of the evolution model 


The tangent linear operator of the curve evolution model around a previous estimate </> required within the 


incremental Gauss-Newton assimilation process reads: 



w — a t A5(j)i + (a n 


/. 


<*t) 


V0fV 2 <tyiV^ V 2 0i 


- 2 - 


V 


V0i 


V4>i 


f. 




V4>i 




The corresponding adjoint can be computed from the adjoint of the discrete scheme used for the tangent 
linear operator. This strategy is straightforward and exact but nevertheless requires storing a lot of intermediate 
information at each grid point when advanced numerical schemes are used. A simpler scheme consists to 
operate a discretization of the adjoint of the tangent linear operator with the same discrete scheme. In our case 
this operator is defined as: 
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(a) (b) (c) 

Figure 2: Differences between two strategies of initialization of the level-sets, (a): a meteosat image 
including a convective cell on the right part; (b) the corresponding initial level-set obtained with a signed 
distance and (c) the corresponding initial level-set obtained based on the temperature. One can see that the 
apparition of a new cell in the left part of the image will be less favored with the signed distance function (b). 
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Both strategies are equivalent when the discrete schemes used to implement the adjoint operator constitutes 
the adjoints of the discrete schemes used for the direct model. This is the strategy we adopted in our imple¬ 
mentation. Exact discrete adjoints operators of the direct model have been used for the adjoint evolution law. 
Some hints on the numerical schemes on which we relied are given in the following. 


4.2.4 Numerical schemes 

In practice, the two level sets of the system state and the velocity field w are defined on the image grid. The 
direct dynamical system (32-33), the backward adjoint equation and the forward incremental model associated 
respectively to operators (39) and (40) are integrated using an explicit discretization scheme with a third order 
Runge-Kutta scheme that respects a Total Variation Diminishing property [44]. The advective term of the level 
set dynamics has been also implemented through a TVD numerical scheme in order to achieve an accurate and 
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stable discretization. We employed a semidiscrete central scheme with a numerical convection flux derived 
from the monotone Lax-Friedrichs flux [25,50]. The second order operators have been implemented through 
traditional centered discrete schemes. 


4.2.5 Observed displacements field 


The dynamical system (39) we consider depends on the local mean of an external displacements field (cf 
relation (38)). This field can be estimated in practice from the image sequence using a fluid motion estimator. 
The estimator on which we relied in this study is formalized as the following minimization problem [12,13,52]: 


arg mm 


/( — 
•' 1 dt 


V/ ■ d + IQ + a(/(|V£|) + /(|V£|)) 


dx, 


(41) 


where I = I(x,t) stands for the image luminance function, d = d(x,t) = (u(x,t),v(x,t)) T is the velocity 
field to estimate, a a smoothing parameter and £ = du/dx + dv/dy (resp. £ = dv/dx — du/dy) is the 
divergence (resp. the vorticity) of the flow. The function /(•) is often the quadratic norm but it can either be a 
“robust” or “softer” penalty function that enables to deal with corrupted data (see [22] for some introduction 
to robust statistics). Roughly, this states that the velocity field to estimate should satisfies 


1. the mass conservation law applied to the image luminance: dl/dt + V/ • d + /£ ~ 0 that assumes 
that the variation on the density is related to the diverging motions. It has been proved in [13] that this 
assumption holds in IR imagery and ; 

2. a spatial consistency that allows the estimation of highly diverging and rotating motion fields, unlike 
usual optical flow techniques based on a first order smoothing of the components u and v. 


Variations of such schemes has been successfully used for the estimation of atmospheric 2D or 3D wind 
fields [13,19,20]. An example of a motion fields estimated through such an estimator on IR images is shown 
on figure 3. 
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Figure 3: Illustration of the estimation of the velocity fields, (a) a convective cell image superimposed with 
the estimated velocity and (b) the velocity field. 

4.3 Initial surface 

The initial surface is in practice often chosen as a signed distance function to the contour. Since a distance 
function is characterized by a unity gradient modulus (||V<^|| = 1), this strategy leads, following the same 
strategy than (36), to a constant covariance at the initial time and hence does not appear to be an efficient 
choice in our case. 

As in our application a simple thresholding technique provides the initial contours, it appears more natural 
to infer the error covariance between the initial state and the initial contour (the background error) from the 
data, and besides to guarantee some consistency with the model error covariance at the initial time. This can 
be simply done fixing the initial surface as a smoothed temperature surface observed at the same moment. 
This smoothing is operated in practice through a Gaussian convolution with a small standard deviation and 
guaranties differentiability properties of this surface. Furthermore, this initialization has the great advantage to 
include in the implicit surface structural informations of the temperature map before an eventual cell formation. 
This avoids in empty cells areas imposing level sets values that depend only on the geometry of the other 
detected active cells - or worst on an artificial boundary when no cell has been detected. Figure 2 shows an 
example of the benefit provided by this initialization. 
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4.4 Observations 


A common way to track objects in images consists in relying on a reference histogram. However in our 
application, such a generic description of the regions of interest is difficult to specify as the spatial organization 
of the temperature inside a convective cell strongly evolves along time. This evolution is depicted in figure 4. 
This figure illustrates a typical evolution of the temperature histograms observed along time for a given of a 
convective cell. From this illustration, it is obvious that a fixed reference object histogram would be strongly 



Temperature (150 K to 250 K) 


Figure 4: Evolution of the non-normalized temperature histogram of a convective cell sequence. Each 
line corresponds to a representation of the temperature distribution at a given instant of the sequence. The 
Y-axis is the successive images of the sequence and the X-axis represents the temperature values observed in 
the cell. 

sensitive to the instant of acquisition of the first image of the sequence. We thus propose an observation 
operator based on the temperature threshold related respectively to the convective cell contour and the cell’s 
core. Assuming, as it is commonly done in convective cells detection issues, that the core’s activity (resp. 
its external area) lies in temperatures under a given threshold T c (resp. above T e ), the function 0, whose 
zero-level corresponds to the border to extract, should be such as: 

H(I(x) - T.)( 1 ^ H (</>{*))) + H(T. - /(x))H(</>(x)) = e, (42) 

where • = {e, c} and II is the Heaviside function (H ( x ) = 0 if x < 0 and H (x) = 1 elsewhere). This relation 
is zero if I(x) > T. and </>(x) > 0 or if I(x) < T. and <j>(x) < 0. This measurement model favors thus the 
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curve interior to be above of a given threshold and the exterior to be lower of the same threshold. Otherwise 
this function is always positive. Hence, this relation enables to define the observation system y = H(^) as: 


J y = o 

[ H(</>) = H(I(x) - T.)(l - H(</>{*))) + H(T. - I{x))H(<f>{x)) 

The operator H is non-linear with respect to (j>. Its associated tangent linear operator yields: 


(43) 




i-H(I(x) - T.) + H(T. - /(*))] m, 


(44) 


where 5(<j>) is the dirac function applied to qb. 


4.5 Error covariance tensors 

The dynamical model error covariance matrices Q and the background error covariance B have been defined 
in section 4.2 and 4.3 respectively. They are account for the quality of both the measured velocities from 
the image sequences and the input luminance data. As for the observation error, a similar procedure has been 
applied: since the observations are extracted from an image threshold independently at each point of the image 
grid, the error observation covariance tensor has been defined as a diagonal matrix whose values are directly 
linked to the surface gradient magnitude (i.e a thresholding measure is associated to high covariances value 
and hence low uncertainties -resp. a low value and high uncertainty- for high magnitude of the temperature 
gradients -resp. low temperature variations-). The proportionality coefficient used in this covariance matrix 
constitutes the single parameter of our assimilation process. In practice this coefficient has been set to 1 in all 
the experiments on which the technique has been run. 

Let us now present some experimental results. 


5 Experimental results 

This section describes some experiments that have been carried out on real data of the Meteosat Second 
Generation satellite. Before showing the results, let us first describe the datasets that have been used as 
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(a) (b) 

Figure 5: Examples of IR 10.8 images on the two different sequences of the benchmark. 

benchmark. 

5.1 Data and geographical area 

Our different image sequences were extracted from August 8th 2005 (time OOhOO) to September 1st (time 
23h59). The upper left hand corner coordinates of the geographical area are (0,-20) degrees in latitude and 
longitude, and (25,40) at the opposite point. The region comprizes a great part of central Africa, which ex¬ 
tends from the Sahel in the northern part of the image to the equatorial region in the southern part. These data 
occurred during the Sahel rainy season, spreading from June to September [2], and during the monsoon rains, 
which takes place in the southern part. The following bands were used: WV 6.2, IR 9.7, IR 10.8, and IR 13.4. 
IR images were transformed into brightness temperatures using equation (1). 

We have selected two sequences, respectively named “single cell” and “multiple cells”. The figure 5 depicts 
one image of two different selected sequences. In the single cell sequence (figure 5 (a)), a single large con¬ 
vective cell is visible. With a diameter of approximately 240 kilometers at its mature stage and a temperature 
below —6067, this convective cell belongs to the category of mesoscale events also called super cells. The 
multiple cells sequence of figure 5 (b) is more complex. It exhibits numerous diagonally organized smaller 
cells, and might be more considered as a multicellular storm than as a true convective cell. Indeed, the biggest 
cell’s diameter does not spread more than one hundred kilometers, and the temperature does not dive much 
below the core temperature threhold. 
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Figure 6: Results on the single cell sequence. The first line corresponds to the contours obtained after the 
propagation of the initial contour with the proposed level set’s dynamics; the second line shows the final results 
after the assimilation process. The contour of the cell is represented in black while the core is in white. 

5.2 Experiments on the single cell sequence 

As explained previously, in order to ensure a non overlapping between the two curves represented by <f> e and 
<f> c (respectively related to the external contour and the core of the cell), the surface <fi c is defined over a 
support corresponding to the interior of the domain delineated by zero level set of <f> e . Its assimilation depends 
therefore on the estimation of surface <f> e . The estimation of the cell boundary surface (be does not depend 
however on the core surface and can be performed first. The assimilation of the core surface is then handled 
afterward. 

The assimilation results are shown in Figure 6. The first line represents five images of the sequence where 
the propagation of the initial contours through the dynamical model are superimposed (the contour in black 
represents the external cell and the one in white is related to the core). The second line depicts the results after 
the assimilation process. On the first row of this figure, one can observe that the dynamical model provides 
already a coarse approximation of the final contour location for the two initial cells (that are merging later in 
the sequence). During the initial integration, we see that the both cells evolve independantly from each other 
along the sequence (first line of Figure 6). The benefit of the assimilation process then clearly appears on 
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the second line of the figure 6 since the two initial cells merge together in a unique cell. It can be observed 
that the estimated curves related to the contour and the core of the cell are well in line with the visualized 
sequence. They also are in accordance with the usual behavior of these convective systems. Let us now study 
the relevance of out approach when several cells split and merge. 

5.3 Experiments on the multiple cell sequence 

This second experiment represents a more complex scene since it contains several cells at different degrees of 
evolution. In addition, as few convective phenomena are initiated in the first image (i.e. only the biggest cells 
are detected through the threshold technique on the first image), the first integration of the dynamical system 
has provided only a rough initial system state for the most proeminent cells: the very first forecast of the initial 
flat temperature surface (first line of Figure 7) corresponding to the initial condition indeed obviously failed 
in predicting a complete series of contours for the different cells observed. The results after the assimilation 
process are shown in Figure 7. Despite an incomplete and unsatisfactory initial solution, we can observe 
that the process enables to estimate accurately the contours and the cores of the different convective cells. The 
different split and merge were correctly managed (the complete sequence can be seen at http://www.sites.univ- 
rennes2.fr/costel/liama/corpetti/ThomasjCorpetti_eng/ANRJ4SG.htm[). 

These two experiments on real meteorological data demonstrate the ability of the proposed assimilation 
process to recover reliable curves describing the contour and the core of the observed convective cell. 

To highlight the robustness of the assimilation procedure, let us now study the relevance of the proposed 
technique in presence of highly corrupted data. 

5.4 Robustness in presence of noisy and corrupted data 

The Meteosat Second Generation data are received, analyzed and stored by the EUMETSAT consortium. For 
some specific applications, it is also possible to install a simple reception station using a conventional TV 
antenna. The MSG signal is captured and transform into a luminance brightness. In such simple reception 
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Figure 7: Results on the multiple cell sequence. The first line corresponds to the contours obtained after 
the propagation of the initial contour with the proposed level set’s dynamics; the second line shows the final 
results after the assimilation process. The contour of the cell is represented in black while the core is in white. 

stations, it is frequent that the signal is corrupted due for instance to interference with aerosols, clouds, ... 
Therefore, many images are likely to be unusable. To analyze the stability of our approach in presence of 
outliers, we have performed the assimilation on a corrupeted version of the single cell sequence where the 
following transformations have been applied: 

• the fourth image has been half bottom cut; 

• a Gaussian-noise with a known standard deviation has been introduced in the seventh; 

• the same noise has been also introduced in the tenth, and the top left-hand quarter has been discarded; 

• the last image has been suppressed. 

On this sequence, the motion field have been estimated using the cost function of relation (41). Thanks to the 
robust norm /(.) used, the outliers areas are “rejected” of the minimization procedure and the solution consists 
in a regularization of the motion field by propagating its divergence and vorticity on the “noisy” data. The 
Figure 8 presents the results after the assimilation superimposed to the original images. The method turns out 
to be very efficient since the estimated curves are still very closed from the ones obtained from the uncorrupted 









Figure 8: Illustration of the assimilation results on corrupted data; the half bottom part of the fourth image has 
been cut out, the seventh has been blurred by a known Gaussian noise, the tenth has a quarter missing and is 
blurred by the same amount of noise, and finally the last image has been suppressed. The estimated level sets 
still correctly delineate the contours of the convective cell and its core, even in the case of missing image 

data. On all the noisy areas, the associated covariance matrices are very high since the observation provides 
no relevant information. As mentioned in previous sections, these latter are indeed related to the magnitude of 
the gradient of the image luminance. Therefore, on such outliers regions, the observation adds no information 
and the system state results in a propagation of the dynamical system. This experience emphasizes well the 
robustness of the presented approach with respect to different potential sources of degradations of MSG data. 

6 Conclusion 

In this paper we have proposed an original and efficient variational assimilation framework for the tracking of 
convective cells from both IR and WV MSG imagery. The technique relies on an optimal control formulation 
which consists to monitor the trajectory of the curves descriptors from the complete sequence of data. Each 
step of the process is composed of two iterative phases: the forward linear development of the state variable 
from a dynamical model, and the backward integration of an adjoint dynamics incorporating a data measure¬ 
ment model. A time consistent family of state variables is by this way generated as a satisfying trade-off with 
respect to uncertainties of the system dynamics and of its observation. Excepted the proportionality coefficient 
used in the background error covariance, all the other parameters involved in the proposed data assimilation 
process are estimated from the image data and do not require any ad-hoc tuning. 
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The cells are characterized using a state variable composed of two levelsets, related respectively to their 
boundary and their core. The level set formulation enables a natural topology adaptability that fits well to the 
possible splitting and merging phases underwent by a convective cell during its life cycle. This approach is 
also relevant to properly process such highly deformable cloud structures. The use of two level sets allows 
identifying the cold clouds from the other surrounding and warmer ones, but also enables to focus on the actual 
convective activity part of the convective cell and to distinguish it from the whole cloud. 

The dynamical model associated to the level sets assumes that the curves are propagated by an external 
velocity field (in practice estimated on the image but that can be provided by any measurement techniques) 
where the uncertainties along the normal and the tangent of the curve are taken into account. This provides a 
generic stochastic process allowing to guide the curves through an implicit surface. Of course, the quality of 
this dynamical model is related to a global accuracy of the velocity fields but local velocity errors are handled 
in our scheme. 

We have outlined the strength of the variational assimilation formulation in its ability to deal with missing 
or altered data. Furthermore, it provides a continuous framework for the state variable tracking, and provides 
therefore a continuous trajectory of each cell. This enables us to provide a cell localization even in case of 
missing or altered images and to oversample the final trajectory (with regard to the image temporal lapse rate). 
In the intricate situation where a cell undergoes multiple merges and splits simultaneously, this decomposi¬ 
tion of the trajectory into shorter intervals of time provides a clearer temporal interpretation of the diverse 
successive topological changes. 

We believe that such a tool is very relevant for the climatologist and meteorology community since the 
analysis of convective cells is still a very intricate problem. 

In order to explore further the potentialities of the level sets assimilation, we present in the appendix a 
tool that labels and estimates several of the geometric and radiative features of the cells that characterizes 
their development’s phase. This is based on several automatic and semi-automatic tools encountered in the 
literature. The association of the assimilation results to the convective cell activity characterization constitutes 
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a valuable and complete tool that might be useful in the meteorologist community. 


Appendix Application to the estimation of characteristic parameters of the 

CONVECTIVE CELLS 

Many studies have demonstrated the interest of characterizing each convective cell individually, as well 
as its phase of development. It is common to characterize the evolution of a convective cell by three states: 
growing, maturation and decay. According to Barrett and Martin [3], each step of the evolution has its own 
particularity. For instance, a cloud produces less rain in its decreasing stage than during its growth. Thus, 
the characterization of the phase should improve the rainfall estimation [32], In addition, several indicators 
related to the stage of development are very useful for specialists. As an example, in [41], the authors have 
presented a technique to estimate precipitations using satellite images based on the “convective cloud area 
evolution’’ index. More generally, several attempts have been carried out in order to provide to the community 
automatic or partially-automatic detection and tracking tools [2,17,18,33]. In this section, we aim at showing 
the interest of the assimilation for the characterization of the life cycle of a given convective cell. From the 
different curves extracted with the proposed assimilation process, we present a tool that allows us to label each 
convective cell, to memorize its history (origin and new cells issued from the current), and to compute several 
significant parameters for the characterization of the convective phenomenon. These parameters are: 

• the averaged brightness temperature; 

• the surface; 

• the perimeter; 

• the averaged divergence; 

• the volume index [46]: J]?ri(/(x) — T) where ni is the number of pixels with the same brightness 
temperature; 

• the gravity center coordinates; 
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Figure 9: Example of labeled convective cells computed from the assimilation results along the sequence. 

• the averaged wind speed vector of the center of gravity, weighted by the temperatures of the cell. 

In addition to the labelisation and characterization of each cell, a “contextual” parameter has been added 
to precise the origin the cell: 

• Normal birth: the cell has just been created and has no antecedent; 

• Normal evolution: the cell has a unique antecedent; 

• Merge: the cell has several antecedents; 

• Split: several cells have the same antecedent; 

These labels allow us memorizing the cell history. For each cell, we are able to identify its date of birth and 
death or its type of evolution (merge, split, ancestor, filiation, ...). To illustrate such labelisation, the figure 9 
represents for the second sequence a picture of the labeled cells at different stages of development. 

The following structure is an example of a cell identity card extracted on this sequence: 

• Label number: 2 

• Born at image #: 4 

• Dead at image #: 7 

• Type of born: NATURAL 

• Futur merge: YES 
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• Futur split: NO 

• Merge with cell(s): 13 

• Overlapping (in pix): 22 

• Average temperature along life: [347.5, 343.5, 340.7, 339.0] 

• Surface along life: [8, 28, 53, 76] 

• Perimeter along life: [8, 24, 35, 39] 

• Average divergence along life: [0.1, 0.2, 0.1, 0.3] 

• Volume index along life: [100,463, 1021, 1593] 

• (Volume index)/(Surface) along life: [12.5, 16.5, 19.3,21.0] 

• Gravity center along life: [(140,46), (140,45), (140,44), (140,42)] 

This expresses that the cell number #2 triggers on image # 4 with a normal birth (it does not result from any 
merge or split of previous cells) and died on image # 7 of the sequence by merging with another cell to create 
the cell # 13. The overlapping between cell # 2 and cell # 13 is 22 pixels. As the cell # 2 is visible on four 
images, four values for each parameter field have been extracted. 

In addition to such an “identity card” for each convective cell, it is also possible to characterize the activity 
of a convective system in the sens of the definition given in [33] that links to the same convective cells future 
merge cells. As already mentioned, we have assumed that a convective cell can be characterized into three 
phases : growing, maturity and decay. Following [39], this convective activity can be characterized considering 
the joint evaluation of the local divergence contained in the 2D motion fields and the temporal variations of 
the brightness temperature. For instance, the growing activity corresponds to the time interval in which the 
averaged temperature of the cell decreases while the averaged divergence is positive; the maturity is associated 
to a calming down activity where the top cloud stops cooling down and its averaged temperature can even 
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(a): bom (b): growing (c): maturity (d): decay (e): death 

Divergence moyenne temperature mean evolution 




Average divergence Average temperature 

Figure 10: Characterization a convective cell’s activity. For a given cell (cell labeled #2 in this figure), an 
image corresponding to each phase and key events is represented. The mean displacement of the cell is also 
superimposed to its gravity center. The second line corresponds to the mean divergence and brightness inside 
the cell. The vertical red lines corresponds to a merge with others cells whereas blue lines represents splits. 

start rising up slightly whereas the decay is characterized by a period where the divergence is negative in 
conjunction with an increasing temperature. Figure 10 illustrates the phase characterization for the cell labeled 
#2. We have also added both indicators (evolution of the mean divergence and temperature inside the cell) that 
enables to characterize its activity. 

With such available data, it is possible to print out other useful characterization indicators. For instance, 
the figure 11 exhibits some of common descriptors used in different previous studies for the characterization 
of convective cells [2,17,18,33]. From these plots, one can see that even if some indicators are correlated 
in some phases (for instance, a positive divergence is related to a positive gradient surface evolution), they 
contain all a specific information. The figure 11 (c) exhibits the quantity — div(i^) where A stands for 
the area of the cell. This index, which links the area evolution to the average divergence has been introduced 
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Figure 11: Area and gradient of the area of the cell 2 of the second sequence along the 23 images of its 
existence. 


in [29], is null during the growing phase. The others indexes, related to the surface, perimeter, volume and 
temperature, are estimated from the estimated curves and the associated motion fields. They may provide 
precious information that can be used for specialists to better analyze and understand the convective cell’s 
behavior. Such summary cannot be robustly extracted without the tracking and denoising process performed 
through the curve assimilation. 
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